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Dynamic critical behavio(u)r of a cluster algorithm for the 
Ashkin-Teller model 

Jesiis Salas^ and Alan D. SokaP'* 

^Department of Physics, New York University, 4 Washington Place, New York, NY 10003 USA 

We study the dynamic critical behavior of a Swendsen- Wang-type algorithm for the Ashkin-Teller model. We 
find that the Li-Sokal bound on the autocorrelation time (Tmt.f > const x Ch) holds along the self-dual curve 
of the symmetric Ashkin-Teller model, but this bound is apparently not sharp. The ratio Ti„t,£/Cjj appears to 
tend to infinity either as a logarithm or as a small power (0.05 ^ p ^ 0.12). 



Monte Carlo simulations near critical points 
are hampered by critical slowing-down: the au- 
tocorrelation time T grows like for a system 
of linear size L at criticality. The traditional lo- 
cal algorithms have a dynamic critical exponent 
z > 2. Cluster algorithms [1,2] can in some cases 
do much better: the Swendsen-Wang (SW) algo- 
rithm for the ferromagnetic Potts model [1] has 
z between and 1, depending on the number of 
states and the dimensionality of the lattice. But 
there is little theoretical understanding of the dy- 
namic critical behavior of SW-type algorithms. 
We have only a rigorous lower bound [3] 



"Hnt.e, Texp > const X Ch 

where Ch is the specific heat, and hence 



■^int,5) -^exp 



> - 



(1) 



(2) 



where a and v are the standard static critical 
exponents. 

Obviously we would like to know whether these 
bounds are sharp: that is, does (2) holds as equal- 
ity or as a strict inequality? Unfortunately, the 
numerical data for two-dimensional Potts models 
are not very conclusive. For the Ising model, the 
bound would be sharp if r grows like a logarithm. 
However, the available data are compatible with 
an exponent z ~ 0.23 [1,4], with a logarithm [5], 
and even with a behavior r ~ L^/^logL [4]. For 
the 3-state Potts model the bound seems to be 
not sharp: z = 0.55 ± 0.03 [3] versus a/v = 2/5. 
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Finally, the 4-state Potts model is very special: a 
naive fit to the data gives z = 0.89 ± 0.05 [3], 
which is smaller than a/v = 1. This anoma- 
lous behavior can be explained if we take into 
account the true leading term of the specific heat 
Ch ~ Llog~^^^ L [6]. This suggests that r might 
have a similar multiplicative logarithmic term, in 
which case the bound would be sharp (modulo a 
possible logarithm). 

There is yet another way of "interpolating" be- 
tween the 2-state (Ising) and 4-state Potts mod- 
els: both are particular cases of the Ashkin-Teller 
(AT) model [7]. As a matter of fact, the self-dual 
curve of the symmetric AT model joins the critical 
points of these two models (see Figure 1). Along 
this curve the static critical exponents vary con- 
tinuously. 

An SW-type algorithm for the AT model was 
first devised by Wiseman and Domany [8]. Here 
we study a simplified ("embedding") version of 
this algorithm: we want to know how the dy- 
namic critical exponent z behaves along the self- 
dual curve, and in particular whether the Li- 
Sokal bound (l)/(2) is sharp. We have performed 
simulations at three different points on the AT 
self-dual curve (see Table 1). In addition, we 
have reanalyzed the data obtained by Baillie and 
Coddington [4,9] for the Ising model at criticality. 
We conclude that the bound (1) is almost but not 
quite sharp: the ratio Tmt.e/C'ij tends to infinity 
as L ^ oo, either as a logarithm or as a small 
power (0.05 <p< 0.12). 

The Ashkin-Teller (AT) model [7] is a general- 



2 



Phase Diagram of 2D AT Model 




Figure 1. Phase diagram of the symmet- 
ric Ashkin- Teller model on the square lattice. 
The self-dual curve is B-DIs(decoupled Ising)- 
P(Potts)-C. The solid curves represent second- 
order phase transitions, and the dotted curve is 
the non-critical part of the self-dual curve. The 
dash-dotted curve is the 4-state Potts-model sub- 
space. The roman numerals designate the differ- 
ent phases of the model (see text). 



ization of the Ising model to a four-state model. 
To each lattice site x we assign two Ising spins 
CTj, = ±1 and Tx = ±1, and they interact through 
the Hamiltonian 



H 



-J (TxCTy- 

(xy) 



- J ^ ] Ta; Ty - 



(3) 



where the sums run over nearest-neighbor pairs 
(xy). Note that the fields cr, r and err play sym- 
metric roles in this model; we can consider any 
two of these three as the "fundamental fields". 
The plane K = corresponds to a pair of de- 
coupled Ising models with interactions J and J', 
while the line J = J' = K is the 4-state Potts 
model with Jpotts = 4J. 

The family (3) of AT Hamiltonians exhibits 
several symmetries. First of all, we can permute 
freely the spin variables (cr, r, err); this is equiv- 
alent to permuting the couplings (J,J',K). Sec- 



Point 


J = J' 


K 


P (4-state Potts) 

ZF 

X2 

DIs (decoupled Ising) 


0.274653 
0.302923 
0.344132 
0.440687 


0.274653 
0.220343 
0.147920 




Table 1 

Points of the self-dual curve of the symmetric AT 
model where our MC simulations were performed. 
We also include the values corresponding to the 
Ising model (DIs). 



ondly, if the lattice is bipartite we can flip cr or 
r or both on one of the two sublattices; this is 
equivalent to flipping the sign on any two of the 
three couplings (J, J', K). 

In this paper we are interested in the 2D 
square-lattice symmetric (J = J') AT model: 

Hs = -j'Y [o-xO-y + TxTy) - (TxTx<TyTy .(4) 

{xy) {xy) 

This model exhibits a rich phase diagram [10,11], 
which is shown in Figure 1. There are four differ- 
ent phases: 

I) Baxter phase. The spins cr and r are indepen- 
dently ferromagnetically ordered. 

II) Paramagnetic phase. Here the three spins cr, 
r and err are disordered. 

III) The spins cr and r are disordered, but their 
product err is ferromagnetically ordered. 

IV) Both cr and r are disordered, but their prod- 
uct err is antiferromagnetically ordered. 

The curves separating these phases are critical. 
All except the curve B-P (i.e. the self-dual curve) 
are expected to be Ising-like, and their exact lo- 
cations are unknown. The self-dual curve is given 

by [11] 

e-^^=sinh2J, (5) 

and the part above the 4-state Potts point P is 
critical. The critical exponents vary continuously 
along the self-dual curve, and their values are ex- 
actly known by relating the AT model with the 
Gaussian model [12]. 

Let us consider the general AT Hamiltonian (3) 
with the condition 



J, J' > \K\ 



(6) 
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(Actually, by using the symmetries, any AT 
model on a bipartite lattice can be mapped onto 
an equivalent model satisfying this condition.) 
One can then construct a SW-type algorithm for 
this model, following the general scheme of [13]: 
the idea is to decompose the Boltzmann weight 
associated with a given bond as a linear combi- 
nation of Kronecker deltas of the spins, and then 
to introduce new auxiliary variables which live 
on the bonds. The final result is essentially the 
same algorithm as introduced by Wiseman and 
Domany [8]. Details of the derivation can be 
found elsewhere [14]. For the "multi-cluster" ver- 
sion of this algorithm, we have proven [14] the 
Li-Sokal bound (l)/(2) by following the scheme 
of the original proof [3] . 

A simpler SW-type algorithm can be intro- 
duced by considering the Boltzmann weight of a 
given bond (xy), conditional on the {r} configu- 
ration: it is 

where p^y = 1 — exp(J + KT^Ty). This system of 
a spins can be simulated using the standard SW 
algorithm with ferromagnetic effective nearest- 
neighbor coupling = J + Kr^Ty. An identical 
argument applies to the {r} spins when the {cr} 
spins are held fixed. The "embedding" algorithm 
for the AT model consists therefore of one SW 
update of the {cr} spins with the {r} spins held 
fixed, followed by SW update of the {r} spins 
with the {cr} spins held fixed. 

We have simulated this embedding algorithm 
for the 4-state Potts model at criticality and for 
two other points (X2 and ZF) on the self-dual 
curve (see Table 1). X2 was already considered in 
[8], and ZF is the model introduced by Zamolod- 
chikov and Fateev [15]. We have studied lattice 
sizes ranging from L = 16 to L = 512 (1024 for 
the 4-state Potts model). The number of mea- 
surements ranges from 8 x 10^ to 4.4 x 10®. (In 
units of the autocorrelation times, our run lengths 
are always at least 10'*r, except for L = 1024 
where we reached only 1500r.) In all cases we 
discarded 10^ iterations for equilibration; this is 
always > 150r. For the Potts case we have mea- 
sured the energy, specific heat, second-moment 
correlation length and susceptibility. In the other 



two cases, we have measured these quantities for 
both the cr, r spins and for the product err. 

We have performed standard weighted least- 
squares fits to a power-law Ansatz for all the di- 
verging static quantities (specific heat and sus- 
ceptibilities). In these fits we use the data with 
L > a cutoff Lmin ) and we vary Lmin until we ob- 
tain a good value; in this way we can protect 
against corrections to scaling. Our results (see 
Table 2) agree fairly well with the exact answers. 
The specific heat of the 4-state Potts model is the 
only exception. But in this case the true leading 
behavior is ~ Llog~^^'^ L [6], not merely ~ L. If 
we include this theoretical input in the fit we get 
a much better estimate [14]. 

We have also measured the autocorrelation 
functions and the corresponding integrated and 
exponential autocorrelation times (see e.g. [16]). 
Power-law fits yield the dynamic critical expo- 
nents Zint,£ reported in Table 2. We see that the 
value of Zint,£ is always slightly higher than the 
effective value of a /v. Thus, the Li-Sokal bound 
(2) is satisfied — but apparently not as equality 
— along the AT self-dual curve (5). 

Nevertheless, the very small values obtained 
here for Zint,£ — ajv (less than 0.11) suggest that 
the bound (l)/(2) might be sharp modulo a log- 
arithm. To check this, we have studied the ratio 
'Hnt.e/C'if- This ratio is in all cases an increasing 
function of L. We tried fits to const + AL~'^ and 
const + j4/logL, but the results are poor: either 
"X^ is too large, or the value of the constant A 
is implausibly large. We conclude that T^^t^^zjCii 
probably does tend to infinity as L ^ oo. The 
question is: in what way? We next tried to fit 
T\w\.,e I Ch to a pure power law AL^ . The values 
are very good in all cases; the power pis small and 
seems to increase from the Ising and X2 models 
(p ^ 0.05) to the 4-state Potts model (p ^ 0.12). 
If this is the true behavior, it would mean that 
the bound (1)/ (2) fails to be sharp by only a small 
power. On the other hand, we have also tried the 
Ansatz Tmt.e/C'ij = A + BlogL. The fits are 
again quite good for all the models. In this case 
the bound would fail to be sharp only by a multi- 
plicative logarithm. It is very hard to distinguish 
numerically between these two scenarios: a loga- 
rithmic behavior can be quite well mimicked by 
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Ratio 


4-state Potts model 
numerical exact 


ZF model 
numerical exact 


X2 model 
numerical exact 


Ising model 
numerical exact 


ll^ 
ol/u 

Zint,£ 


1.744 ±0.001 7/4 
1.744 ±0.001 7/4 
0.768 ± 0.009 1 X log-^/^ 
0.876 ±0.012 > 1 X log-^/^ 


1.750 ±0.004 7/4 
1.668 ± 0.005 5/3 
0.663 ±0.006 2/3 
0.740 ± 0.010 > 2/3 


1.751 ±0.001 7/4 
1.605 ±0.001 1.6045 
0.438 ±0.008 0.4183 
0.477 ±0.028 > 0.4183 


7/4 
1/2 
log 

0.240 ± 0.004 > log 



Table 2 

Static critical exponents and dynamic critical exponent (zint.e) coming from power-law fits to the Monte 
Carlo results [14]. For the Ising model we include our own fits to the dynamical data reported in Refs. [4,9]. 
Errors are one standard deviation. The symbol 1 x log~^^^ (resp. log) means that the leading term of the 
specific heat for the 4-state Potts model (resp. the Ising model) behaves like Llog~^^^L (resp. logL). 



a power law when the range of variation of log L 
is not very large. Indeed, we can equally well fit 
the data to a function log^ L. Details of the data 
and the fits will be reported in [14]. 

In summary, the fits of the ratio Tmt.e/C'ij sug- 
gest two different scenarios 

,„ \ ALP with 0.05 <p< 0.12 

To distinguish between these two behaviors would 
require high-precision data on significantly larger 
lattices than those simulated here, probably up 
to at least L = 2048. 

These computations were carried out at NYU, 
the Pittsburgh Supercomputing Center, and 
the Cornell Theory Center. The authors' re- 
search was supported in part by a M.E.C. 
(Spain)/Fulbright fellowship (J.S.), and by NSF 
grant DMS-9200719 (J.S. and A.D.S.). 
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